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Abstract 

Streamer electrical discharges are often investigated with computer simulations 
of density models (also called drift-diffusion-reaction models). We review these 
models, detailing their physical foundations, their range of validity and the most 
relevant numerical algorithms employed in solving them. We focus particularly 
on schemes of adaptative refinement, used to resolve the multiple length scales in 
a streamer discharge without a high computational cost. We then report recent re- 
sults from these models, emphasizing developments that go beyond cylindrically 
symmetrical streamers propagating in homogeneous media. 



1. Introduction 

Streamers are transient electrical discharges that appear when a non- 

conducting medium is suddenly exposed to a high electric potential. While the 
average background field might be too low for plasma generation by electron im- 
pact on neutral molecules, the streamer discharge channel can enhance the electric 
field at its growing tip so strongly that it can create additional plasma and prop- 
agate nevertheless. The streamer achieves this through a very nonlinear dynam- 
ics with an intricate inner structure and locally very steep density gradients. This 
structure generates the local field enhancement and maintains propagation. There- 
fore, the accurate simulation of single, cylindrically symmetric streamer channels 
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is far from trivial, even if the electrons are approximated by their densities as is 
conventionally done. 

It should be noted that the full process evolves on even more scales. On the 
one hand, when the electric field ahead of the streamer head is very high, single 
electrons can run away and generate hard radiation through Bremsstrahlung. To 
describe these effects, single electrons have to be followed in a particle model; 
numerical methods to track these particles in an efficient manner and the derivation 
of density models are discussed in this special issue by Li et al. [|5|]. On the 
other hand, in the laboratory, in technological applications and in lightning related 
processes, hundreds to ten-thousands of streamer channels can propagate next to 
each other. In this case it is vital to develop models on a more macroscopic scale 
than the density approximation. Through model reduction, one can derive moving 
boundary models llq-llOll for the underlying ionization fronts, or one even can try 
to develop models for the streamer channel as a whole. 

Within the current short review, we discuss the justification for density models 
for streamers, the numerical solution strategy, and then some results that go be- 
yond single cylindrically symmetric streamers in homogeneous media. We treat 
interacting streamers, first through the trick of considering an infinite array of 
identical streamers next to each other, then through full 3D simulations. Next, 
we discuss streamer branching in full 3D. Finally we discuss how to simulate the 
emergence of sprite streamers from earth's ionosphere; a peculiar issue is here 
how to combine the very different scales of the electric fields generated by the 
thundercloud-ground-ionosphere system with the fine inner structure of the dis- 
charge and the varying density of the atmosphere. 

With the density models, we here focus on the oldest and most extensively 
investigated family of numerical streamer models; they are also known as fluid 
models, continuous models or reaction-drift-diffusion models. These models are 
now mature enough that their predictions often deviate from experimental obser- 
vation by a margin comparable to the random variations and uncertainties of the 
experiments. 



2. Model formulation 

2.1. The density model 

A classical density model for a streamer discharge has the structure of a reaction- 
drift-diffusion equation for the electrons and reaction equations for various ions 
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and excited species coupled to the electric field 

= V-(n ef i e E) + V-D e Vn e +S? + Sf, (la) 

= sr+sf\ i=i,...,N, (lb) 

£ V • E = q, E = -V0. (lc) 

Here n e is the electron number density and [Z,] is the density of the heavy species 
Z, denoting an ion or an excited state, and p, e and D e are mobility and diffusion 
coefficient of the electrons. S'") denotes mostly local generation or loss of species 
due to reaction at direct encounter of particles; the most prominent example is 
electron impact ionization, i.e., the generation of electron-ion pairs through im- 
pact of electrons on neutrals; the efficiency of this process strongly depends on 
the electron energy that in turn is determined by gas density and electric field. 
S^l denotes the generation of species through radiative transport in a generically 
nonlocal process; the most prominent example in streamer physics is the gener- 
ation of electron-0 2 pairs through photo-ionization in air. Finally, E is the elec- 
tric field, (p the electric potential, and q is the local space charge, determined as 
1 = [Z,] - e n e ; here e is the elementary charge and q, = -e, 0, -l-e is the 

charge of species Z,. 

2.2. The underlying approximations 

We list here the main approximations underl ying this model. For a more 
through derivation of the model, we refer to Ji, 11, 12]. 

1 . The electrons are much more mobile than the ions, and the dynamics takes 
place on their time scale, therefore ion motion is neglected — whether the evalu- 
ation time scales exceed this limit, needs to be checked after evaluations. 

2. The electron motion is approximated by drift and diffusion within the local 
field. This entails that the electrons rapidly relax to a velocity where the acceler- 
ation by the electric field exactly cancels the momentum losses due to collisions 
with other particles, and that spatial or temporal variations of the electric field 
are not important. The accuracy of the drift-diffusion approximation and the local 
fiel d ap proximation for electron currents in streamer ionization fronts was verified 



in 111211 for electric fields up to some threshold; this statement holds for arbitrary 
gas density. However, when the field at the streamer head exceeds a limit (180 
kV/cm in air at standard temperature and pressure, STP |0, [HI [HI, [3]), some 
electrons in the tail of the electron energy distribution function might not relax to 
a steady velocity anymore, but run away with a field dependent probability. Above 
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the thermal runaway electric field (260 kV/cm in STP air 11151-11711) the bulk of the 
electrons will keep accelerating up to relativistic energies and the drift approxi- 
mation breaks down. 

3. The field is calculated in electrostatic approximation. For typical current 
densities and diameters of streamers Rial , this approximation is justified. 

4. The ensembles of electrons and ions can be treated in density approxi- 
mation. The validity of this approximation depends first on the gas density, and 
second on the specific region within the streamer. First, streamers in different gas 
densities n are related through similarity relations or Townsend scaling in very 
good approximation 11181 . Il9ll ; these relations imply that the electron density n e in 
a streamer scales as n 2 , and that the intrinsic length scales scale as 1 jn. Conse- 
quently, the total number of electrons in similar parts of streamers scales as 1 /n. 
This implies that the density approximation in all regions of the streamer becomes 
better when the gas density decreases. On the contrary, the density approximation 
for electrons and ions becomes worse when the density of the neutral medium 
increases, and it has to be reconsidered at the high densities of liquids. (Other 
approximations to be reconsidered in streamers at liquid densities are the absence 
of heating and the neglect of electron-electron and electron-ion collisions in the 
source terms 5„„; in gases only electron-neutral collisions are included in typi- 
cal streamer models, but the degree of ionization and therefore the probability of 
electron-electron or electron-ion collisions increases with increasing density be- 
cause n e /n oc n.) A review of validity and corrections to these similarity solutions 
is given in II 1811 . 

Second, different spatial regions of the streamer have to be distinguished: in- 
terior, front and outer region. The streamer interior has large densities of electrons 
and ions and little inner structure; the density approximation is certainly justified 
there. In the outer region, electron densities are typically very low (except for 
high discharge repetition frequencies or with other strong ionization sources) and 
densities should rather be interpreted as probabilities; however, a change of this 
low background ionization has essentially no influence on streamer propagation in 
typical situations Q20Q , but it can have an influence on branching and "feather" for- 
mation in very pure gases with little nonlocal photo-ionization II20L 12111. Finally, 
the streamer ionization front connects interior and exterior and the electron den- 
sity can develop steep gradients, in particular, in high electric fields and with little 



trp 

or no nonlocal photo-ionization lU2ll20l-l22ll. Within this region, the local density 



approximation can break down as electrons with a higher energy are ahead of oth- 
ers. The electron energies are then determined not only by the local electric field, 
but also by the gradient of the electron density; this effect can be incorporated in 
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a modification of the impact ionization term; or alternatively, the model can be 
extended by an electron energy equation [[Till . These terms can be included in the 
general structure of the reaction-drift-diffusion model above and do not require 
any separate numerical strategy. According to 111 in . the density approximation 
describes streamer propagation well, as long as no massive electron runaway sets 
in. However, to properly describe the effect of density fluctuations and electron 
run-away on streamer stability, a particle approach should be chosen. Presently 



available models are Monte Carlo particle models rtl6ll23l . 12411 and spatially hybrid 



models |5, 11] that combine particle and density models in different volumes of 



the simulation. 

2.3. Typical models in air 

The general structure of the density model defined above has been recognized 
decades ago but the search for better reaction and transport parameters in different 
gases goes on. Most modeling work is devoted to air modeled as a mixture of 
nitrogen and oxygen. 

The impact ionization term S " n in air includes the generation of free electrons 
and positive ions through electron impact on neutral molecules as well as elec- 
tron losses due to attachment to oxygen. In particular the first reaction strongly 
depends on the electric field. There is also field dependent detachment of elec- 
trons from fl25|], however, it is typically neglected in air, and it indeed only 



contributes in a relevant manner in nitrogen with admixtures of oxygen as low as 



1 p.p.m., but not at the 20% oxygen content in air [20J]. 



This is because air has very strong nonlocal photo-ionization — for a his 



torical review of this concept, we refer to the introduction of [21]. In oxygen 



nitrogen mixtures such as air, photo-ionization is in most cases included through 



the Zhelezniak model Il26ll . In this model, the states h 1 !^, b' 1 !^ and c^D* of 
molecular nitrogen, populated by impact excitation in the high-field regions of 
the streamer, decay by emitting photons, some with energies above the ionization 
threshold of molecular oxygen, 12. 1 eV (1025 A). The Zhelezniak model assumes 
that the excitation of the emitting states is roughly proportional to the impact ion- 
ization. Then the production of electron-ion pairs per unit of volume and time is 
written as 

c * p « f /<Plr-r1)Sf(rV(pr') 

Sph(r) = - . : -r 7 , (2) 

4np + p q J \pr-pr'\ 2 

where £ is a proportionality constant, p q = 60 Torr = 80 mbar is the quenching 
pressure of the emitting states and h is the absorption function of photo-ionizing 
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radiation. Photo-ionization might also be very strong in hydrogen-helium mix- 
tures that dominate the atmospheres of gas giants like Jupiter and Saturn 112711 : 
however, quantitative models are not available for this mixture. 

The reaction rates strongly depend on electron energies (that are typically pa- 
rameterized by the electric field that accelerates the electrons), but also the trans- 
port coefficients /u e and D e depend on it, though weakly. In most literature, they are 
treated as constant. Furthermore, in general the diffusion coefficient D e is a ten- 
sor with different components parallel and perpendicular to the electric field 
but most investigators approximate it as a scalar. Reaction and transport coeffi- 
cients are mostly calculated in advance from cross-section data such as ll2~8ll by 
solvin g th e Boltzmann equation for electrons with numerical codes such as BOL- 
SIG+ Ir29ll . Alternatively, they also can be calculated by performing Monte Carlo 
simulations with good statistics IH, U]. Much work is currently devoted to im- 
provements of the cross-section data as well as to the evaluation of the Boltzmann 
equation |3()| . 

On the other hand, it can be stated that streamer discharges seem to be amaz- 
ingly robust both against model changes in simulations I120n as well as against 
changes of gas composition in experiments 112 ll 12711. The reason is probably the 
strongly nonlinear dynamics of the streamer that approaches generic attractors of 
the dynamics independently of microscopic details. E.g., in streamers in nitrogen 
with 1 p.p.m. oxygen, the effect of background ionization or photo-ionization can 
hardly be distinguished, and major changes of the photo-ionization model in air 
have only minor effects on streamer velocity [|20ll . 



3. Numerical schemes 

Generally speaking, the physical models that underlie numerical streamer sim- 
ulations are always small variations on the model discussed above. These varia- 
tions include e.g. the use of field-dependent mobility and diffusion and the inclu- 
sion of kinetic equations for additional chemical reactions. 

To reduce computational demand, almost all streamer simulations are per- 
formed in cylindrical symmetry with radial and longitudinal coordinates (r, z) al- 
though below we will discuss some recent extensions to real 3d. 

3.1. Discretization of the transport equation 

The main challenge in the numerical solution of the densities of charged par- 
ticles in streamers arises from the very steep gradients in the streamer tip. The 
choice of discretizations has been largely driven by the requirement of handling 



6 



these gradients efficiently. The first studies 113 ll - l34ll used variations on the Flux- 
Corrected Transport (FCT) algorithm developed by Zalesak [35], combined with 
a modified Euler scheme for the time integration. In 1995 Kulikovsky ll36ll in- 
troduced a scheme based on improvements of the Scharfetter-Gummel algorithm 
that has subsequently been used also by Liu and Pasko B7I1 and Bourdon and 
coworkers 113 811 . On the other hand, first-order upwind schemes were used in Il39ll 
but they seemed to be too diffusive, and they over- stabilized the streamer fronts. 
J 2211 . an upwind-biased scheme was applied that used the Koren flux limiter 
to switch between first-order upwind and third order upwind-biased. 
These discretizations may be combined with various algorithms to work with 
different resolutions in different volumes of the domain. The importance of this 
part of the numerical scheme stems from the multi- scale nature of the streamer 
process discussed above. Kulikovsky 114 ill used a "window" with a fixed, hi gher 
resolution that followed the streamer head, while Pancheshnyi and coauthors ll42ll 
and Eichwald and coauthors I143ll used a rectangular product grid with smoothly 
varying space steps. Min and coauthors I144ll and Nikandrov and coauthors [45] 
have used dynamically adaptive meshes. 

Montijn and coauthors [|22ll used a scheme of nested uniform grids generated 
adaptively with a criterium based on the second derivatives of the species and 
charge densities as follows. The refinement builds a hierarchy of grids starting 
from the coarsest grid that covers the entire volume and descending to finer grids. 
To decide the areas to refine, at each grid at level /, monitor functions M u are 
calculated, where u stands for the electron density or space charge density. These 
functions come from the discretization of 



M u (r,z) 



(Ar 1 ) 2 



1 d I du 
r dr\ dr 



+ (Az 1 ) 2 



d 2 u 



dz 2 



(3) 



where Ar 1 and Az' are, respectively, the /-level grid resolutions in r and z. The re- 
finement criterium is then specified from grid-independent refinement tolerances 
e u as 

refine all grid cells i, j where — > e u . (4) 

max w . 

u 

Streamer fronts are pulled, i.e. the dynamics is determined in the leading edge of 
the ionization front Il46ll . For curved fronts as occuring in streamers, the leading 
edge is limited by the decay of the electric field. Therefore the refinement cri- 
terium © is often insufficient in the leading edge of the front. The approach used 
in I122ll is to use a region larger than given by © and extend all grids in the leading 
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edge up to the layer where the electron density is below a small threshold. How- 
ever, when photo-ionization is included in the simulations the electron density 
ahead of the front decays much slower and this approach leads to an overexten- 
sion of the refinement grids. Instead one can use a criterium based on the absolute 
value of the electric field, adding to © the additional criterium 

refine grid cells i, j where \E t j\ > E threshM , (5) 

with an appropiately choosen E thresho i d . 

3.2. Strategies for the Poisson-equation 

Solving the Poisson equation is usually one of the most expensive parts of 
simulations of density models for streamers. It is challenging, because the solution 
is nonlocal and instantaneous, i.e., the boundary conditions together with the very 
localized charge distribution in the simulation domain determine the electric field 
instantaneously. The Poisson equation is most often discretized in a cartesian grid 
by a five-point second order scheme. The resulting linear system is then solved 
using either iterative methods, usually from the Succesive Over-Relaxation (SOR) 



family II37L 143. 147L 14811 or by direct methods such as cyclic reduction 11221 14911 
and SuperLU pQll . In general, direct methods can be much faster than iterative 
methods for these problems, but they perform similarly if the iterative method is 
provided with a good initial guess from the solution of a previous timestep. 

Even with fast algorithms, solving the Poisson equation is a major computa- 
tional bottleneck for streamer simulations in full three dimensions. Two earlier 3D 



streamer simulations 11251. bill constitute a proof of principle that such computa- 



tions can be performed. An efficient and parallelizable approach for that problem 



is described in 115211 : it extends the grid refinement technique from 112211 to 3D sim- 
ulations. This method is efficient for streamer dynamics that can be represented 
in cylindrical coordinates (r, z, 6) with a relatively low resolution in the axial co- 
ordinate AO = 2n/N. 

The discrete Fourier transform (DFT) in 6 of the electrostatic potential <p(r, z, 9) 
reads 

N-l 

Mr, z) = Yj <Kr, z, 2nnlN)e- 2 « iknlN ', (6) 

n=0 

which implies the usual properties of the DFT of a real quantity with N even, 

<p k (r, z) = ffir-k( r > Mr, z) = fair, z), <p N /2(r, z) = <ftf /2 (r, z), (7) 
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where * indicates complex conjugation. The benefit of working in Fourier space 
is that different modes k are decoupled and can be solved in parallel. The Laplace 
operator V 2 can be decomposed into V 2 = V 2 Z + p- -§52 an d the midpoint discretiza- 
tion of the second derivative with respect to 6 is 

cb } , -2<A, + <*,-., 
D wi ~ ~^0i ' W 

where for clarity we have omitted the indices for r and z. Applying the DFT as 
defined in © we obtain 

D g (p k = — 2 = — Iwjtl (p k , (9) 



where 



\w k \ 2 = -^(l-coskA9). (10) 



[V 2 J0,(r, z) - l ^ffa(r, z) = -q k {r, z)/e , (11) 



Note that as AO — > 0, \wk\ 2 — > k 2 , giving the well known form of the continuous 
Fourier transform of a second derivative. 

Finally, the equation for including the r, z terms is a two-dimensional Helmholtz 
equation 

r* 

where [V 2 .] is the discrete Laplace operator in (r, z) and qk(r, z) is the DFT of the 
space charge density. The available numerical methods for this equation are basi- 
cally the same as for the Poisson equation (where the second term is missing) but 
(fTTI) makes it possible to solve the equations independently and hence in parallel 
for each k. 

Unfortunately, the reaction terms in (fTal)-([Tbl) are nonlinear so the convection- 
diffusion-reaction side of the model must be solved in real space, where it can 
also be trivially parallelized. Hence an efficient parallel algorithm must switch 
between Fourier and real space at each timestep. This is relatively fast by using 
widely available Fast Fourier Transform (FFT) codes. 

Since the Poisson equation is usually solved in a cartesian grid, a second prob- 
lem concerns the modeling of a curvilinear or point-shaped electrode. This is an 
important aspect of streamer simulations since in experiments streamers almost 
always emerge from a pointed electrode. One common approach, introduced by 
Babaeva and Naidis I153n . is to consider a spherical needle electrode where the 
computational domain of the densities begins only at the tip of the electrode, and 
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to calculate its influence by the method of image charges. Another approach was 
proposed in 154(1 . based on a simplified version of the method of virtual charges. A 
similar approach that does include the needle electrode into the domain of density 
computations is the Charge Simulation Method II55L 15611 that was implemented for 
streamers in section 7.2 of [1570. Recently the use of immersed boundaries was 
introduced in streamer simulations by Celestin et al. [50y. 

A non-planar electrode can be also simulated by u sing curvilinear coordinates, 
adapted to the geometrical shape of the electrodes 115 8l. 15911: this is particularly 
appropriate if the pointed-electrode has a parabolic shape. 



3.3. Strategies for the nonlocal photo -ionization 

In the first simulations of streamers, photo-ionization was not included in the 
model but rather substituted by an unrealistically high constant pre-ionization [31]. 
Later, the Zhelezniak model was introduced by directly evaluating the integral 
©, i.e., by summing the effect of each emitting source on each grid point in 
the simulation domain. The number of computations needed in this approach 
scales quadratically with the total number of grid points, growing much faster than 
the number of computations required for the Poisson equation or the advection- 
diffusion time stepping. One can use a coarser grid for photo-ionization calcula- 
tions 114 ill , but this seems to introduce spurious effects. Another approach, used in 
[60], is to approximate the spatial distribution of the photo-ionizing emitters by a 
sphere centered at the center of mass of radiation and with a radius equal to some 
characteristic radius of the streamer. 

A more accurate approach is to approximate the Zhelezniak integral by the 
solution of elliptic partial differential equations. This method was introduced 



independently by Segur et al. 116 ill and Luque et al. [62J], and further refined 
iniMSE]. 



4. Beyond the cylindrically symmetrical streamer in a homogeneous gas 

The keystone of streamer simulations consists of a single cylindrically sym- 
metrical streamer propagating in a homogeneous medium. As we mentioned 
above, this is by itself a physically relevant and challenging numerical problem. 
Nevertheless, in recent years several extensions to this paradigm have been intro- 
duced that investigate a wider range of problems. Here we review three of them: 
the interaction between neighboring streamers, the branching of a streamer chan- 
nel and the propagation of streamers in the upper atmosphere, with a significant 
gradient of air density. 
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4. 1 . Streamer interaction 

Streamers very rarely appear in isolation; most often they belong to bunches 
of many streamers, either because they have emerged at different points of a wire 
electrode or because they are contained in a highly branched tree emerging from a 
point electrode. Since each streamer carries charge, they interact electrostatically. 

4.1.1. Bunches of streamers and Saffman-Taylor streamers 

The first work to address the interaction of neighbouring streamers was by 



Naidis 116511 . He used a quasi-2D density model to estimate the effect of the 
interaction from neighbouring streamers in a regular one-dimensional lattice of 
well-separated streamers. The underlying assumption was that the radius of each 
streamer is much smaller than the distance to the closest neighbour. 

A further step was demonstrated in [66], where the transversal dynamics of 
interacting streamers is taken into account by a modification of a single-streamer 
model. As shown in Fig. Q] (left) the interaction between neighbouring streamers 
is incorporated as a lateral homogeneous Neumann boundary condition of the 



electrostatic potential. As a first step, Ref. 16611 was limited to a 2D (planar) model 



of negative streamers without photoionization. However the technique can also be 



applied to more realistic situations H67I1 . 

The outcome, after a transient regime where each streamer expands, is a steady 
state where each channel fills half of the available width L (the distance between 
one streamer and its neighbour). In this steady state, the electric field far behind 
the array of streamer heads is completely screened. The numerical solution shows 
a remarkable agreement with the Saffman-Taylor analytical solution of the well- 



known problem of viscous fingering in Hele-Shaw cells (Fig.[H right) II68L 16911 . 

The analogy between streamers and viscous fingering is explicit in moving 
boundary models 0-[l(l H • 



4.1.2. Interaction between two streamer heads 

Another approach to the problem of streamer interaction was presented in I152ll . 
where the method of 3d calculations explained above was used to study the inter- 
action between two negative streamer heads propagating in parallel in a homoge- 
neous electric field (see Fig.0. 

In general, for a given separation between the streamer heads, the outcome of 
the interaction depends on the gas composition. In cases where photo-ionization 
is absent such as in pure nitrogen (as far as it can be realized in the laboratory 



112 111 ), the electrostatic repulsion between the charges in the streamer heads drives 
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Figure 1: Left: periodic array of negative streamers. The interaction between neighboring stream- 
ers can be incorporated into a density model by imposing homogeneous Neumann boundary con- 
ditions on the electrostatic potential on the dashed lines. Right: The resulting streamer profile fits 
the analytical solution of a selected SafFman-Taylor finger, with width L/2. (reproduced from [66]) 




Figure 2: Left: Scheme of refinement for the 3d simulation approach described in the text. The 
red surfaces surrounds the streamer head, defined by its charge density. The slices represent the 
numerical resolution of the density equations: they are equi-spaced in the angular direction but in 
the r, z plane they contain a scheme of nested grids, common for all slices. Right: Surfaces of 
constant electron density of two interacting negative streamers in nitrogen at atmospheric pressure 
propagating downwards in a homogeneous electric field of 80 kV/cm. Reproduced from 15211 . 
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them appart, bending the streamer channels outwards. However, for gas compo- 
sitions where photo-ionization plays an important role, such as air, the ionization 
in the volume between the two heads counteracts the electrostatic repulsion be- 
tween them, driving the two streamers to merge. In general the outcome of the 
streamer interaction would therefore depend on the separation between them and 
the balance of these two counteracting processes. 

4.2. Branching 

The branching of a streamer channel is perhaps the most challenging phe- 
nomenon for numerical simulations. The main reason is that branching is caused 
by a dynamical instability at the streamer front, where small irregularities in the 
very thin space-charge layer grow and eventually form new streamers. Therefore 
simulations of the branching process demand a high numerical accuracy to rule 
out numerical artifacts. This difficulty is compounded by the fact that branch- 
ing completely breaks the cylindrical symmetry of a single streamer channel, thus 
requiring full 3D simulations for the complete understanding of the phenomenon. 



In 2002 Arrayas et al. 117 111 studied the branching of a negative streamer in a 
minimal model that did not include photo-ionization. There, cylindrical symme- 
try was imposed and the result was interpreted as an upper bound for branching. 
The resulting picture was one of branching due to a Laplacian instability, in close 
analogy to other well studied physical models such as diffusion-limited aggre- 
gation and viscous fingering. Following this result, the numerical scheme was 
improved by Montijn et al. [22], who showed that branchingtime converged as 



the highest resolution of the numerical grid was improved 1172(1 . thus ruling out 
that branching was due to numerical artifacts. 

Although those results imposed a cylindrical symmetry on the streamer, we 
can show now that they are robust when that restriction is removed. We use the 



method of 3D simulations detailed in 115211 to investigate the dynamics of a neg 



ative streamer without photo-ionization under small perturbations that break its 
cylindrical symmetry. We study here only a linear regime in which the solution of 
our system of equations remains at all times close to a cylindrically symmetrical 
solution; in a realistic situation the deviations from this solution are likely larger 
and the streamer will sooner deviate significantly from the symmetrical solution. 

In our simulation a negative streamer propagates in nitrogen (i.e. there is no 
photo-ionization nor electron losses via attachment) between parallel plates in a 
14.4 mm gap with a homogeneous field of 50 kV/cm. For the spatial discretization 
we use the nested-grids approach of Montijn et al. ll22ll . discussed above. The 
highest resolution in r and z is Ar = Az = 2.5 /jm 
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First we run a simulation with perfectly symmetrical initial conditions con- 
sisting in a semispherical gaussian neutral ionization seed attached to the upper 
electrode. The radius of the seed is 73.6 /urn and its highest electron and ion den- 
sity is 2.6 • 10 10 cm~ 3 . As seen in Fig. [31 the streamer branches between 6.75 
and 7.5 ns. In this context, "branching" means that the streamer front loses its 
convexity. 

To show that this branching is indeed associated with the breaking of cylin- 
drical symmetry observed in experiments, we run a similar simulation where the 
initial condition deviates slightly from cylindrical symmetry. In this simulation 
we used cylindrical coordinates with angular resolution AO = In IN, N = 32. 

The asymmetrical initial state is constructed from the symmetrical initial elec- 
tron densities n e (r,z) by drawing 771, ... , 77^-1 numbers from a normal distribution 
N(0, 1) and then setting 



where n ek (r, z, t) is the discrete Fourier transform of the asymmetrical electron 
density, defined as in ©: 



The strength of the perturbation to the symmetrical state is e and the factor r/R 
is introduced to keep n e (r, z, 9) continuous at r = 0. The typical length R is 
introduced to keep e dimensionless. Here we take Ro = 2.3 fim (the inverse of the 
Townsend multiplication rate), e = 10~ 8 but we also used e = 10 6 to check that 
we were indeed in a linear regime and the growth rates of the perturbation do not 
depend on e. We chose to perturb only the initial electron density, keeping the 
initial ion density symmetrical. 

Note that this way of perturbing the electron density does not guarantee that it 
will remain positive. However, here we are interested in very small perturbations; 
in this limit the resulting density will almost certainly be positive everywhere and 
we checked for that. 

To track the small deviations from perfect symmetry observed in the asym- 
metrical simulation we define 



n e0 (r, z,t = 0) = n e (r, z) 

r 

n ek {r, z,t = 0) = erj k n e (r, z)— (k t 0), 



(12) 
(13) 
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Figure 3: Evolution of an axi-symmetrical streamer head represented by contours delimiting the 
region where charge density is larger than half of its maximum. The contours are plotted at regular 
timesteps of 0.75 ns. Note that the equivalent figure for the evolution of a streamer with a slight 
perturbation of the initial symmetry looks almost identical to this one due to the smallness of the 
perturbation. 
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where u stands for either the electron density n e or the charge density n, - n e . The 
W[ u) give us an idea of the "spectral content" of electrons or charges in each mode. 
Note that Wy' e) is the total number of electrons and W^ 1 '""^ is the total charge 
contained in the simulated volume. In general the W ( k u) are. N complex numbers 
k = 0, . . . , N - 1 but since they come from the DFT of a real quantity, they satisfy 
relations analogous to ©; for N even they reduce to (N - 2)/2 complex numbers 
and 2 reals. To analyze the spectral content of the electron density we look at the 

amplitudes \Wf\ = yJw { k u) W ( k u) * = ^W ( k u) wf_ k for k = 0, ... , N/2. 

Still, plotting |W£ (01 is usually not the best way to visualize the evolution 
of the streamer. First, the number of electrons and the total charge changes in 
time; one should compare each mode with the corresponding total, \Wq \. Second, 
because we give random initial values to each \Wi |, it is also better to normalize 
them according to their initial value. Hence we define 

Vl"\t) = * ° - , (16) 

K\t)\\w?W 

which satisfies V ( Q u \t) = 1 for all times, and V ( k \0) = 1 for all k. The evolution of 
the V k '\t) is shown in Fig.gl 

The evolution of the axial perturbation clearly shows a change of regime at 
t = t B with t B ~ 7 ns, coinciding with the onset of streamer branching. For t < t B 
the evolution of the Fourier modes is smooth and all of them decay with respect to 
the k = mode (dV k ' e) /dt < 0); during this phase the streamer is becoming more 
symmetrical. On the other hand for t > t B some modes k > 1 grow with respect 
to the k = mode. This is also shown in Fig. [51 where we represented the growth 
rates y of the unnormalized components W? 1 ' of the space charge density; for 
each k, y is obtained by fitting Wr'(t) to an exponential Ae 7 ' for t < t B and t > t B . 

These results show that "branching" interpreted as non-convexity of an axi- 
symmetric solution is indeed related to faster growth of the axial non-symmetric 
modes relative to the k = mode and hence to "branching" in the more conven- 
tional sense of bifurcation of a single channel. 

We note that there are two possible reasons for this outcome: (a) that the 
streamer reaches a state of general dynamical instability that simultaneously af- 
fects both the cylindrically symmetrical and the asymmetrical modes and (b) that 
a non-convex, cylindrically symmetric streamer head is unstable against perturba- 
tions of its symmetry; this is expected if the electric field has a significant off-axis 
peak. These two processes are compatible and likely both play a role. 
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Figure 4: Time dependence of the normalized Fourier components of the electron density [VT (i); 
see text for a definition] (left) and the space charge density [Vr 1 : (f)] (right) in a simulation 
with a slightly asymmetrical initial condition. Note that the V) (t) with k > decay smoothly 
up to t — ts ~ 7 ns. This means that the k — component of the electron density (i.e. the total 
number of electrons) is growing faster than any other Fourier component; we interpret this as the 
streamer becoming closer to cylindrical symmetry. At t — fg a change in behavior is clear; after 
this points several v5 u '(t) with k > grow faster than the k — mode and the cylindrical symmetry 
is gradually broken. 
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Figure 5: Growth rate of the unnormalized Fourier components of the space charge density 
[ffj ], We plot separately the growth rates before and after the branching time tg ~ 7 ns. 
Note that for t > tg several Fourier modes grow faster than the Icq mode, indicating that the shape 
of the streamer deviates from cylindrical symmetry. 
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4.3. Streamers in inhomogeneous media: sprite streamers 

Up to now, we have only considered streamers propagating in homogeneous 
media. Streamers models have been extended to inhomogeneous gas densities 
and compositions in order to study the interaction between a streamer and gas 
bubbles [73] and the propagation of a streamer inside a plasma jet IT74n . Here we 
consider another important case where streamers propagate in an inhomogeneous 
gas, namely sprites [75Q. These are transient discharges in the upper atmosphere 
above active thunderclouds, where air density varies with altitude. Sprites span 
altitudes from about 50 km up to 90 km, are tens of kilometers wide and usually 
last for some milliseconds. Often they are composed of two separate regions: an 
upper, diffuse one usually extending from about 80 to 90 km altitude and a lower, 
filamentary region composed by hundreds of streamers IU9L 17611 (sometimes called 
"sprite tendrils"). Pasko has recently reviewed models of single sprite streamers 
propagating in homogeneous air density lf77ll . 

A single streamer in the filamentary region can extend for tens of kilometers 
in altitude. In this range the density of air varies considerably (roughly by a factor 
of 2 every 5 km). Typical streamer lengths approximately scale with density IU8I1 
so the density variation introduces another length scale on top of those discussed 
above. The first numerical models to include the variation of air density in sprite 
streamers were published recently 11781. 17911 . 

The geometry of sprite simulations is depicted schematically in Fig. |6j The 
electrical field is created when a thunderstroke at an altitude L « 10 km moves 
a large charge to the ground. Most sprites are generated after positive cloud-to- 
ground lightning strokes; these strokes leave a negative charge Q in the cloud. In 
the quasi-electrostatic approximation fl8Qil Q is taken to change so slowly in time 
that the electromagnetic radiative terms can be neglected (i.e. the speed of light 
is infinite). The upper and lower boundaries of the resulting electrostatic prob- 
lem are, respectively, the Earth's surface and a horizontal layer in the ionosphere 
at about 85-90 km altitude, where the Maxwell relaxation time due to the atmo- 
spheric conductivity is much shorter than the characteristic times involved in the 
discharge. 

In principle, the charge Q can be included in the simulation domain and treated 
as a prescribed source term in the Poisson equation. However, the distance be- 
tween this charge and the upper boundary of the simulation is about 80 km and in 
order to minimize the effects of the artificial lateral boundaries the domain would 
have to be extended to a radius much larger than these 80 km. To avoid this 
overextension, in 117 811 the field created by Q was calculated semi-analytically as 
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-90 km 




Figure 6: Scheme of the geometry of sprite streamer simulations. First the field of the cloud 
charge with boundary conditions on ground and ionosphere is solved analytically. The the Poisson 
equation for the sprite charges is solved in the complete domain between ground and ionosphere; 
the diffusion-drift-reaction equations are solved only in the upper, shaded volume. 
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the field created by a charge at a distance L from the lower of two parallel elec- 
trodes. This is expressed as a sum of the fields created by an infinite number of 
image charges; this sum is truncated after 4 terms (2 on each side). To this one 
adds the self-consistent field of the sprite, calculated by solving the Poisson equa- 
tion with homogeneous Dirichlet (tp = 0) boundary conditions in the upper and 
lower electrodes. The horizontal width of the domain must be larger than the lat- 
eral extension of the space charges to minimize the effects of artificial boundaries 
but it can still be smaller that the distance between Q and the sprite. 

Once the field is calculated, the diffusion-drift-reaction equations must only 
be solved in a smaller volume (see Fig©. In 117 811 the volume was a cylinder from 
55 to 85 km in altitude, with a radius of 20 km. 

A notable difference between streamers in the laboratory and sprite streamers 
is that in a laboratory they are almost always initiated from a sharp metallic elec- 
trode that is clearly not present in the upper atmosphere. Raizer 117 611 proposed that 
streamers emerge from "plasma patches" created by electro-magnetic waves 118 lh . 
Recently, Marshall et al. 118211 used a finite difference, full electromagnetic code 
to evaluate the influence of lightning-generated electromagnetic pulses (EMP) on 
the electron density in the upper atmosphere, reporting large enhancements of up 
to a factor 4 for repeated EMP's. 

However, an isolated patch would create a two-headed streamer propagating 
simultaneously up and down [13711 . In observations, however, sprite streamers al- 
ways start propagating downwards 11831 . 18411 . In 117 8ll it was shown how a sin- 
gle, downwards-propagating sprite streamer emerges from the destabilization of 
a sharp screening-ionization wave from the ionosphere; figure [7] shows the opti- 
cal emissions resulting from that process. The initial condition amounted here to 
air density exponentially decreasing with altitude and electron density increasing 
with altitude as in fair weather night-time. 

A relevant recent result related to sprite streamers concerns the trailing optical 
emissions emanating from the streamer body at a certain distance from the head. 
Although initially this luminosity was attributed to chemical processes in the wake 
of the streamer head [851 lS6|], two recent papers 11791 18711 used streamer density 
models to show that it originates from a re-enhancement of the electric field. This 
ocurs if there is a substantial electrical current in the streamer body, generated 
because the net charge in the streamer head increases during propagation. In ll79ll 
the electric field in this re-enhancement approached E « E^, this appears to be 
dynamically stable in a current-carrying state. 
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Figure 7: Optical emissions emitted from a streamer emerging from the sharpened edge of a 
relaxation-ionization wave from the ionosphere. This wave produces diffuse optical emissions 
observed as a halo. (Reproduced from [78]) 



5. Outlook 

As with any scientific endeavour, it is hard to predict the future of streamer 
density models. The ultimate aim of simulating and understanding a fully branched 
tree of streamers strictly within the model described in this paper looks unrealistic 
in the short and medium term. 

On the other hand, gradual improvements to density models will surely come 
in the short run. There are two obvious directions to extend the physical mod- 
els. In the "microscopic direction" the modeling of the microscopic phenomena 
is improved; this includes e.g. a more accurate calculation of the photo-ionization 
sources, additional chemical processes and the inclusion of stochasticity, converg- 
ing towards detailed particle and hybrid models. The "macroscopic direction," 
extends the simulation domain, possibly to include many interacting streamers 
and also incorporate in an improved fashion the interaction between the streamer 
and its surroundings. Although in principle they are compatible, both directions 
demand a higher computational it is unlikely that they will be combined in a prac- 
ticable unified model. Even if approached separately, they will likely require faster 
computers and numerical algorithms. 

However, it may not be necessary to describe the full range of physics, from 
the microscopic to the macroscopic level, inside a single model. A more promis- 
ing approach is to develop a hierarchy of models of decreasing microscopic detail. 
Each model would rely on simplifications and assumptions extracted from models 
down in the hierarchy or from direct experimental observations. Density models 
play a pivotal role in this hierarchy, sitting between particle models and higher- 
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level descriptions of a streamer tree such as the dielectric breakdown model of 
Niemeyer et al. ll88l l89ll. phenomenological branching models [90] and moving 
boundary models if^l— Tl~oh . 
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